Sobre os dados

Focos de queimada

Os dados de focos foram obtidos do banco de dados do INPE: https://queimadas.dgi.inpe.br/queimadas/bdqueimadas#exportar-dados

Os dados foram filtrados pelos limites Amazônia legal, entre os anos de 2011 e 2020.

Outras variáveis

Tanto a distância à rodovia mais próxima quanto a distância à área urbana mais proxima foram calculadas em função de shapefiles de rodovias e áreas urbanas disponibilizados pelo IBGE. Os dados de áreas desmatadas são elaborados pelo INPE através do PRODES.

Imports

Os pacotes sf e raster são essenciais para a manipulação de dados vetoriais e matriciais respectivamente. O pacote tidyverse contém funções uteis para a manipulação de tabelas. O pacote spatstat é utilizado para extrair estatisticas espaciais como a densidade de focos. Por ultimo, o pacote viridis é usado só para definir as cores nas imagens e gráficos.

library(raster) 
library(tidyverse)
library(sf)
library(spatstat)
library(viridis)
library(randomForest)

Criação da camada focos/km²

A ideia geral aqui é criar uma imagem raster com a contagem da ocorrência de focos por km². Primeiro os dados de foco devem ser carregados. Eu aproveito para ler o shapefile da amazônia legal também.

focos <- st_read("focos_2011_2020.shp")
## Reading layer `focos_2011_2020' from data source 
##   `F:\Pablo\Documentos\Doutorado\side\focos_2011_2020.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1215884 features and 15 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2828897 ymin: 8004499 xmax: 6111664 ymax: 10573720
## Projected CRS: SIRGAS 2000 / Brazil Polyconic
am_legal <- st_read("amazonia_legal.shp")
## Reading layer `amazonia_legal' from data source 
##   `F:\Pablo\Documentos\Doutorado\side\amazonia_legal.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2794537 ymin: 8004219 xmax: 6112012 ymax: 10586380
## Projected CRS: SIRGAS 2000 / Brazil Polyconic

Em seguida criar um raster vazio de resolução de 10km para preencher os dados:

focos_km2 <- raster(ext=extent(focos), crs=crs(focos), resolution=10000)

E enfim contar os focos por pixel da imagem:

focos_km2 <- rasterize(focos, focos_km2, fun='count', field=1, background=0) %>%
  mask(am_legal)
plot(main="Focos brutos por 10km²", focos_km2, col=viridis(256))

Ou então fazer uma análise de densidade:

pts <- st_coordinates(focos) %>%
  unique()
pts_ppp <- ppp(pts[,1], pts[,2], window=as.owin(am_legal), check=FALSE)

A densidade é diferente da quantidade bruta de focos por pixel no sentido de que ela é suavizada em função de uma banda (distância). A função padrão de suavização é a gaussiana.

Um teste rapido usando diferentes bandas:

bandas <- c(10000, 25000, 50000, 100000)
dens <- lapply(bandas, function(x) raster(density(pts_ppp, sigma=x))*1e6)
dens <- stack(dens)
names(dens) <- c("10km", "25km", "50km", "100km")
plot(dens, col=viridis(256))

Nesse caso uma banda de 50km pareceu apropriada.

pt_dens <- density(pts_ppp, sigma=50000, eps=10000)
pt_dens <- raster(pt_dens)*1e6
crs(pt_dens) <- crs(focos_km2)
pt_dens <- resample(pt_dens, focos_km2)
plot(main="Densidade de focos por km²", pt_dens, col=viridis(256))

Lado a lado:

r_stack <- stack(focos_km2, pt_dens)
names(r_stack) <- c("Focos por 10km²", "Densidade de focos")
plot(r_stack, col=viridis(256))

É importante lembrar que os valores calculados representam um histórico de 10 anos, portanto o valor anual seria a média, ou o valor apresentado dividido por 10.

Finalmente, salvar os resultados:

writeRaster(focos_km2, "focos_10km.tif", options=c("COMPRESS=DEFLATE", "PREDICTOR=2"), overwrite=T)
writeRaster(pt_dens, "densidade.tif", options=c("COMPRESS=DEFLATE", "PREDICTOR=2"), overwrite=T)

Variáveis independentes

Os mapas dessas variáveis ja foi feito usando o QGIS, então é só carregar eles.

dist_rod <- raster("dist_rod.tif") %>%
  mask(am_legal)
plot(main="Distância à rodovia mais próxima (km)", dist_rod/1000, col=viridis(256), maxpixels=20000)

dist_desm <- raster("dist_desm.tif") %>%
  mask(am_legal)
plot(main="Distância à área desmatada mais próxima (km)", dist_desm/1000, col=viridis(256), maxpixels=20000)

Amostragem

A princípio a amostragem vai ser feita por via de pontos aleatórios.

set.seed(1)
amostras <- st_sample(am_legal, 5000) %>%
  st_sf('ID' = seq(length(.)), 'geometry' = .)
ggplot()+
  geom_sf(data=am_legal, col="gray", fill=NA)+
  geom_sf(data=amostras, col="red", pch= 3, alpha=0.3)+
  theme_bw()

Extração dos dados rasteriais e conversão para tabela:

dados <- tibble(
  focos_10km2=raster::extract(focos_km2, amostras),
  densidade=raster::extract(pt_dens, amostras),
  dist_desm=raster::extract(dist_desm, amostras),
  dist_rod=raster::extract(dist_rod, amostras)
) %>%
  drop_na(densidade)
dados

Traduzindo em forma de gráfico:

ggplot(dados, aes(dist_rod/1000, densidade))+
  geom_point(alpha=0.3)+
  ylab("Focos por km²")+
  xlab("Distância à rodovia (km)")

ggplot(dados, aes(dist_desm/1000, densidade))+
  geom_point(alpha=0.3)+
  ylab("Focos por km²")+
  xlab("Distância ao desmatamento (km)")

Regressão

Para começar, um teste usando RandomForest para fazer uma regressão:

modelo_rf <- dados %>% 
  randomForest(densidade ~ dist_rod + dist_desm, data=.)
modelo_rf
## 
## Call:
##  randomForest(formula = densidade ~ dist_rod + dist_desm, data = .) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##           Mean of squared residuals: 0.02922064
##                     % Var explained: 41.8
dados2 <- dados %>%
  mutate(preds = predict(modelo_rf, newdata=.), erro=(densidade-preds))

Comparação com a realidade:

dados2 %>%
  ggplot(aes(dist_rod/1000, densidade))+
    geom_point(aes(color="Observado"),alpha=0.3)+
    geom_point(aes(y=preds, color="Previsto"), alpha=0.3)+
    xlab("Distância à rodovia")+
    ylab("Densidade")

dados2 %>%
  ggplot(aes(dist_desm/1000, densidade))+
    geom_point(aes(color="Observado"),alpha=0.3)+
    geom_point(aes(y=preds, color="Previsto"), alpha=0.3)+
    xlab("Distância ao desmatamento")+
    ylab("Densidade")

Distribuição de erros:

dados2 %>%
  ggplot(aes(x=erro))+
  geom_histogram(bins=50)+
  scale_x_continuous(breaks=seq(-1,1,0.1))

Enfim classificar a imagem:

dados_stack <- stack(dist_desm, resample(dist_rod, dist_desm)) %>%
  brick()
pred_rast <- predict(dados_stack, modelo_rf)
plot(pred_rast, col=viridis(256))

par(mfrow=c(1,2))
plot(main="Previsto", pred_rast, col=viridis(256))
plot(main="Observado", pt_dens, col=viridis(256))

Alguns pontos onde é possível melhorar: